Artificial square ice and related dipolar nanoarrays 
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We study a frustrated dipolar array recently manufactured lithographically by Wang et al. [Nature 
439, 303 (2006)] in order to realize the square ice model in an artificial structure. We discuss models 
for thermodynamics and dynamics of this system. We show that an ice regime can be stabilized 
by small changes in the array geometry; a different magnetic state, kagome ice, can similarly be 
constructed. At low temperatures, the square ice regime is terminated by a thermodynamic ordering 
transition, which can be chosen to be ferro- or antiferromagnetic. We show that the arrays do not 
fully equilibrate experimentally, and identify a likely dynamical bottleneck. 



Introduction: The ability to manipulate constituent 
degrees of freedom of condensed matter systems and their 
interactions is fundamental to attempts to advance our 
understanding of the variety of phenomena presented to 
us by nature. For a long time this has been achieved by 
utilizing the combinatorial richness of the periodic table 
of elements to construct different chemical compounds. 

A more recent option is to use the tools of nanotechnol- 
ogy to custom-tailor degrees of freedom which can be as- 
sembled in a highly controlled manner; e.g. this has been 
proposed for realizing a topologically protected quantum 
computer using Josephson Junction arrays |2]. Submi- 
cron superconducting rings have also been used to pro- 
vide effective spin- 1/2 degrees of freedom |2|. 

Very recently, Wang and collaborators have used litho- 
graphic techniques to create a periodic two-dimensional 
array of single-domain submicron ferromagnetic islands 
[l|, depicted in Fig.Q] This design approach takes advan- 
tage of well-established lithographic techniques and en- 
ables reading of the state of the system with local probes, 
such as magnetic force microscopes, to image the state 
of single constituent magnetic islands 4] . 

The first aim of this study is to assemble a system that 
realizes the square ice model. This is an attractive target 
model because of its long and distinguished history dur- 
ing which algebraic correlations and a finite entropy at 
zero temperature have been established, as well as con- 
nections to exact solutions, quantum magnetism, unusual 
dynamics and gauge theories 5]. 

The pioneering study by Wang et al. raises a num- 
ber of important questions which we try to address here. 
Firstly, what are appropriate models for the arrays' ther- 
modynamics and dynamics? Secondly, what other sys- 
tems can one hope to build with these techniques? And 
thirdly, what are interesting directions in which further 
developments would be desirable? 

In particular, the question of whether a dipolar sys- 
tem with long-range interactions can be modelled by the 
short-range ice model is rather similar to the one posed 
in the case of (three-dimensional) dipolar spin ice, where 
it was found that a nearest-neighbor description was sur- 
prisingly accurate for a range of properties, such as the 



FIG. 1: (color online) Left: Atomic force microscope image of 
an array studied in Ref.Jj. The islands have length / = 220nm, 
width 80nm and thickness 25nm. Right: Map of the ratio 
J2/ Ji of the second to the first nearest neighbor interactions 
(highlighted in the left part) for different values of lattice 
constant, a, and sublattice height offset, h. In the white zone, 
I J2/J1 — 1| < 5%. In the left (blue) region, the ordered state 
is antiferromagnetic, whereas it is ferromagnetic in the right 
(yellow-red) area. 



low-temperature (Pauling) entropy j6l. 171. I81I91I1C 

In this paper we show that an analogous equiva- 
lence between ice states and the ground states of two- 
dimensional dipoles on the links of the square lattice is 
more delicate. However, it can be established via a route 
quite different from the three-dimensional case, namely 
by (a) placing the dipoles pointing in different directions 
onto slightly different heights and (b) manufacturing the 
dipoles as elongated as possible. It turns out to be easier 
to realize kagome ice in a dipolar array, as the requisite 
symmetry is compatible with embedding in a plane. 

As a byproduct, the low-temperature antiferromag- 
netic (in ice language: antiferroelectric) instability of the 
original model can be designed to be replaced by a ferro- 
magnetic one. However, the experiments observe no or- 
dering transition, implying that the array does not fully 
equilibrate. We are thus led to study a phenomenologi- 
cal model for its dynamics: zero-temperature ('greedy') 
stochastic dynamics subject to an energy barrier for spin 
flips. This reproduces the experimental measurements 
semi-quantitatively. Such dynamics are insufficient to an- 
neal out isolated defects violating the ice rules even for 



strong interactions, thus preventing the estabhshment of 
an ice — or indeed an ordered — configuration. 

In the remainder of this paper, we first analyze — by a 
mean-field theory and using Monte- Carlo simulations — 
the equilibrium statistical mechanics of arrays designed 
to mimic square and kagome ice. We then turn our atten- 
tion to the dynamics of the arrays studied in Ref. 1 and 
present our results for the local correlations. We conclude 
with a brief discussion of disorder and an outlook. 
Dipoles on the links of a square lattice: The in- 
teractions between the magnetic islands are dipolar, and 
therefore essentially a geometric property of the array 
(lattice constant: a), described by the Hamiltonian 
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where Vij = Vi — r^, and fl{ri) is the dipole moment at 
n, which points along the link. We treat the dipoles 
either as points {l/a -^ 0) or as uniform, monodomain 
thin needles of finite length /. 

Can such an arrangement be used to access square ice 
physics? In other words, is there a (low-temperature) 
regime for such a dipolar system, where configurations 
obeying the ice rule are overwhelmingly present with ap- 
proximately equal weights? 

This is, by construction, the case for the Tsing-ice 
model', in which the four islands emanating from a given 
vertex of the square lattice interact equally and antifer- 
romagnetically, so that in the ground state two dipoles 
point into each vertex, and two out. These are precisely 
the ice states. The Fourier transform of this interaction 
has two branches (as there are two sites in the unit cell, 
one for each link direction of the square lattice), which 
have eigenvalues 

Ml) = ; Uq) = 2J (sin2 I + sin^ |) . (2) 

It is the flatness of the lower branch that indicates the 
frustration of the ice model. As all the ice states are 
linear combinations of the flat-band eigenvectors only, a 
suflicient condition for Eq.^to yield an ice regime is that 
its lower branch be flat and share the eigenvectors of the 
Ising-ice model [l^, at least approximately. 

By contrast, for a nearest-neighbor interaction only, 
one obtains a symmetric pair of bands: 



JiAq) = T2J 



sm — sm -^ 
2 2 



(3) 



The spectrum for artiflcial ice (Fig. ^ lies in between, 
with a substantial dispersion to the flat 'ice' band, of 
about a third of the total bandwidth. Thus, there is 
a thermodynamic transition to antiferromagnetic order, 
which in our Monte- Carlo simulations occurs at temper- 
ature Taf = 1.70(5) Ji (where Jn is the strength of the 
n*^-neighbor interaction, see examples in Fig. ^. 





FIG. 2: (color online) Left: Spectrum for artificial ice (cut off 
after 10a, for point-like dipoles). Right: same for l/a = 0.7 
and h/a = 0.207: the lower band becomes almost flat, (less 
than 1.5% of the total bandwidth). 



The principal reason for the dispersion of the lower 
band is the inequivalence — as in the F-model — of the six 
vertices of the ice model, which fall into two groups. One 
pair (labelled Type I in Ref. [l|, see Fig. ^ has zero total 
magnetic moment, while the others (Type II) have a net 
moment along a diagonal and are higher in energy. This 
inequivalence results from the fact that, unlike the case 
of a tetrahedron in (i = 3, the six bonds between the four 
islands belonging to a vertex are not all equivalent. 

However, this can be remedied by introducing a height 
displacement h between magnetic islands pointing in the 
X- and ^/-directions. 

The ratio J2/J1 of the two inequivalent bond energies 
(Eq. nj is shown in Fig. Q] as a function of h/a and l/a. 
There is a set of choices for these parameters such that 
the interaction energies are approximately equal. For 
point-like dipoles {l/a -^ 0), this value is 



/^ice/a = V (3/8)'/' -1/2 ^0.419, 



(4) 



and taking into account the flnite extension of the dipoles 
lowers the required height offset. For instance, for l/a = 
0.7, hice/ci ^ 0.207. In principle, for 1 — l/a = e ^ 0, 
hice/d ~ V^e -^ 0. However, in this ideal limit the effects 
of disorder, flnite transverse width and a possible internal 
structure of the dipoles will all play a role. 

Having flxed the short-distance trouble by introducing 
the modulation in height, the question remains what hap- 
pens to the long-distance part of the dipolar interaction, 
which in d = 3, amazingly, turned out to leave the ice 
regime intact [3,|9|. However, the mechanism responsi- 
ble for this equivalence in d = 3 10] is not operational in 
d = 2, as it requires the dimensionality of the dipolar in- 
teraction to coincide with that of the underlying lattice. 
Here, however, we have a d = 3, 1/r^ dipolar interaction 
in a d = 2 array. Nonetheless, the present situation is rel- 
atively benign, as the Fourier sum of a 1/r^ interaction 
in d = 2 is absolutely convergent (obviating the need for 
an Ewald sum). 

Further neighbor terms can be suppressed parametri- 
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FIG. 3: (color online) Ising-ice vs. dipolar arrays: frequency of 
vertex types (% deviation from random distribution for single 
vertex, which is |, ^, \ and | for Types I-IV, respectively), 
entropy (S), and heat capacity (C) from Monte-Carlo simula- 
tions; x-axes were scaled for high-T asymptotics to coincide. 
The fraction of non- ice- rule vertices is below 1% for T < 0.42 J 
for all depicted systems. At low T, the ice regime, which 
widens with increasing //a, is terminated by ferro- (dashed) 
or antiferromagnetic (dots) order for different h. 



cally in the ideal limit of l/a -^ 1. The ratio Jn>?>/ Ji,2 
vanishes as e ^ 0, thus yielding the ideal Ising-ice model. 
As l/a is reduced, the flat band initially acquires only a 
small dispersion. To demonstrate this, in Fig. [2 we have 
plotted the mode spectrum for l/a = 0.7, which corre- 
sponds to the a = 320nm sample ^ . The overlap of its 
eigenvectors with those of the Ising-ice model differs from 
1 by less than 0.1% over the entire Brillouin zone. 

This demonstrates that an ice regime can be obtained 
by this route. Our Monte-Carlo simulations on the dipo- 
lar and the Ising-ice model for l/a = 0.7 (Fig. 13 bear out 
this statement: the intermediate ice regime is terminated 
at high T by thermally activated defects violating the ice 
rules, and at low T by an ordering transition. Choosing h 
on either side of the optimal value /iice, this transition is 
ferro/antiferromagnetic, respectively, and perhaps even 
to a more complicated state very close to hice- 
Kagome ice: The ground states of antiferromagnetic 
Ising spins on the kagome lattice define what is known as 
kagome ice, with the 'ice rules' requiring each triangle to 
have one spin pointing in and two out, or vice versa 11]. 
Since the sites of the kagome lattice correspond to the 
links of the honeycomb lattice, one can pose the question 
whether a dipolar array forming a honeycomb lattice will 
display a kagome ice regime. 

The case of kagome ice has the advantage that the 
three bonds of the triangle are equivalent, unlike the 
six bonds of the square. This means that the nearest- 
neighbor Hamiltonian does not require any fine-tuning 
through a height offset. Furthermore, using the above 
limit of e ^ 0, we can again parametrically suppress the 
importance of further- neighbor interactions, and hence 
obtain a representation of kagome ice. The vestiges of 
the further-neighbor terms will again give rise to an or- 
dering transition, terminating the ice phase on its low-T 
side. We note that kagome ice is a phase distinct from 
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FIG. 4: (color online) Frequency of vertex types (as in Fig.l^. 
Interaction energies scale approximately as 1/a^, a factor of 
30 between the extremal points. Experiments (symbols) are 
shown against dynamics simulations for needle dipoles (Eq.Q] 
dotted), and for increased J2 (dashed, see text). 



square ice in that its long- wavelength theory is different — 
its correlations are not algebraic, but exponentially short- 
ranged even at zero temperature. 

Dynamics and annealing: Given the impossibility of 
thermally equilibrating the array, the authors of Ref. |l| 
used a rotating magnetic field Be^t^ gradually stepped 
down, to speed up the dynamics 12|. The question 
whether such an 'algorithm' can be efficiently used to 
find the ground state of a system has been discussed in 
the context of spin glasses 13j. (However, in the case of 
an ice regime, the question is somewhat simpler, namely 
whether it is possible to find one of exponentially many 
ice configurations.) 

For our model of dipolar needles, the energy scales for 
the arrays studied in Ref. 1 are: 

'Zeeman' energy : |/i • 5ext| < 2.6 x lO^K (5) 

'Exchange' : 3.6 x lO^K < Ji < 1.1 x lO^K 

Antiferromagnetic long-range order should thus be 
present at room temperature. 

We consider a phenomenological model for the dynam- 
ics, motivated by the experimental protocol. Firstly, with 
the experimental temperatures well below the interac- 
tion strengths, we use a form of zero-temperature Monte- 
Carlo dynamics. Secondly, we note that vertices violating 
the ice rules are experimentally present throughout. As 
such defect vertices can be removed using single spin flips 
only, we impose a constraint on our single-spin-fiip dy- 
namics: for a flip to be accepted the energy gain must be 
above a threshold 0. Thirdly, the rotating field is mod- 
eled by a field of random orientation. 

This model has two free parameters: the threshold 
and the speed with which the field is ramped down Hi. We 
fix these parameters to be the same for all arrays, and fit 
them to obtain the best agreement with the experimental 
measurements of the local correlations (longer-range cor- 
relations in the experiment are very weak) [l|. The best 
fit is obtained for 6> = 2.3 x 10^ Kelvin and /^ = 87 Tesla"^ 
attempted flips per spin during ramp-down (Fig.EJ. 



This algorithm gives semi-quantitative agreement with 
experiment over a range of interaction energies differing 
by a factor of 30 (Eq. 0). However, we systematically 
overestimate the frequency of Type I vertices compared 
to those of Type II (as do other algorithms we have stud- 
ied). This appears to be due to our needle model (Eq.0 
overestimating the ratio J1/J2. Denoting Sij = %^Z.%\ ? 
where Ei is the energy of an isolated Type i vertex, Ref. 1 
finds from a finite-element simulation: £^2 > 2, £^42 > 6 
for l/a = 0.7. If we reduce the value of | Ji — J2I from Eq.Q 
by 30%, our values grow from £^32 = 1.7, £^42 = 4.9, to 
2.25 and 7, respectively. The resulting fit (Fig.^ dashed 
lines) appears noise-limited |l4J . 

Even for the strongest interactions, about 25% (non- 
ice-rule) defects. Type III vertices, persist. Whereas it is 
easy to remove pairs of appropriately oriented neighbor- 
ing defects by flipping the spin which joins them, such 
an annihilation process occurs with low probability once 
these defects are sparse: they first need to diffuse around 
until they encounter a partner. 

Disorder: Our dynamical model does not take into ac- 
count disorder, which is expected to have substantial in- 
fluence on the dynamical behavior even of single islands 
y . Thus, the good agreement of our model with the ex- 
periment might be due to its correct reproduction of the 
dynamical bottleneck, and not the detailed microscopic 
dynamics: the inability of the defects to 'find' one an- 
other may e.g. simply be due to their becoming pinned. 

Disorder also impacts the ice regime thermodynam- 
ically, as the size of the leading perturbation sets the 
scale for its termination at low temperature. Especially 
for fine-tuned h ^ /iice, disorder might dominate (over 
Ji 7^ J2 and further- neighbor interactions), by selecting 
some ice configurations over others; strong disorder might 
even lead to the presence of defects at any temperature. 

To push the analysis further in this direction, experi- 
mental input would be desirable. What is the variance of 
the islands' geometrical properties? Are there some is- 
lands that freeze at much higher fields than others? Are 
defects usually located at the same positions, and what 
is their spatial distribution? How do correlations evolve 
during the ramp-down of the external field? 
Summary and outlook: We have presented models for 
the dynamics and thermodynamics of frustrated dipolar 
arrays, including ways of stabilizing ice regimes. Perhaps 
the most interesting direction of further study involves 
their dynamics, in the presence of varyingdegrees of dis- 
order. In particular, can other protocols [jj, e.g. involv- 
ing the use of AC magnetic fields, be used to speed up the 
dynamics? Note that in present samples, the largest in- 
teraction energies are more than two orders of magnitude 
above room temperature, so that further miniaturization 
is possible without resorting to cryogenics. 

Better equilibration might then open the door for an 
experimental study of constrained classical 15] and per- 
haps eventually even quantum dynamics 3] (and quan- 



tum ice ll^)- Even though this will require a substantial 
experimental effort, there appears to be no fundamental 
obstacle to obtaining at least a classical ice regime. 

The results obtained along the way should provide in- 
sights into how the physics of frustration can lead to 
new ways of effectively suppressing interactions between 
neighboring magnetic nanoislands and into the limits im- 
posed by disorder, a topic of interest with view to appli- 
cations in memory storage \^. We are thus optimistic 
that dipolar nanoarrays will provide an interesting field 
for further studies. 
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